plot.topo <- function(x,
                      water.z,
                      water.col,
                      water.lty,
                      water.lwd,
                      land.z,
                      land.col,
                      land.lty,
                      land.lwd,
                      legend.loc="topright",
                      asp,
                      mgp=getOption("oce.mgp"),
                      mar=c(mgp[1],mgp[1],0.5,0.5),
                      debug=FALSE,
                      ...)
{
    if (!inherits(x, "topo")) stop("method is only for topo objects")
    opar <- par(no.readonly = TRUE)
#    on.exit(par(opar))
    par(mgp=mgp, mar=mar)
    dots <- list(...)
    if (missing(asp)) {
        if ("ylim" %in% names(dots))
            asp <- 1 / cos(mean(range(dots$ylim, na.rm=TRUE)) * pi / 180) # dy/dx
        else
            asp <- 1 / cos(mean(range(x$data$latitude,na.rm=TRUE)) * pi / 180) # dy/dx
    }
    zr <- range(x$data$z, na.rm=TRUE)

    ## auto-scale based on data in window, if window provided
    if (!is.null(dots$xlim) && !is.null(dots$ylim)) {
        xf <- (dots$xlim[1] <= x$data$longitude) & (x$data$longitude <= dots$xlim[2])
        yf <- (dots$ylim[1] <= x$data$latitude) & (x$data$latitude <= dots$ylim[2])
        zr <- range(x$data$z[xf, yf], na.rm=TRUE)
    } else {
        zr <- range(x$data$z, na.rm=TRUE)
    }

    ## kludge to prevent whitespace above/below or to the left/right
    ## of plots, and also to prevent going past the poles
    lat.range <- range(x$data$latitude, na.rm=TRUE)
    lon.range <- range(x$data$longitude, na.rm=TRUE)
    plot(lon.range, lat.range, asp=asp, xaxs="i", yaxs="i", type="n", xlab="", ylab="", axes=FALSE, ...)

    ## Kludge to prevent latitudes beyond poles
    usr <- par("usr")
    if (usr[3] < -90) usr[3] <- -90
    if (usr[4] > 90) usr[4] <- 90
    if (usr[3] < lat.range[1]) usr[3] <- lat.range[1]
    if (usr[4] > lat.range[2]) usr[4] <- lat.range[2]
    if (usr[1] < lon.range[1]) usr[1] <- lon.range[1]
    if (usr[2] > lon.range[2]) usr[2] <- lon.range[2]

    lines(usr[1:2], rep(usr[3],2))
    lines(usr[1:2], rep(usr[4],2))
    lines(rep(usr[1],2), usr[3:4])
    lines(rep(usr[2],2), usr[3:4])

    xlab <- pretty(usr[1:2])
    if (debug) cat("xlab=",xlab,"\n")
    if (debug) cat("lon.range=",lon.range,"\n")
    xlab[xlab > lon.range[2]] <- NA
    xlab[xlab < lon.range[1]] <- NA
    if (debug) cat("xlab=",xlab,"after trimming\n")
    axis(1, at=xlab, pos=usr[3])
    axis(3, at=xlab, pos=usr[4], labels=FALSE)

    ylab <- pretty(usr[3:4])
    ylab[abs(ylab) > 90] <- NA
    ylab[ylab > lat.range[2]] <- NA
    ylab[ylab < lat.range[1]] <- NA
    axis(2, at=ylab, pos=usr[1])
    axis(4, at=ylab, pos=usr[2], labels=FALSE)

    contour(x$data$longitude, x$data$latitude, x$data$z,
            levels=0, drawlabels=FALSE, add=TRUE,
            col="black")                # coastline is always black
    legend <- lwd <- lty <- col <- NULL
    if (zr[1] < 0) {
        if (missing(water.z)) {
            if (zr[2] > 0) {
                water.z <- pretty(c(zr[1], 0))
                water.z <- water.z[water.z!=0]
                #cat("water.z=");print(water.z)
                ## Do some tricks to get shelf water as well as deep
                if (max(water.z) == -1000)
                    water.z <- c(water.z, -500, -250, -100, -50)
                else if (max(water.z) == -500)
                    water.z <- c(water.z, -400, -300, -200, -150, -100, -50)
                #cat("after tricks, water.z=");print(water.z)
            } else {
                water.z <- pretty(zr)
            }
            water.z <- rev(sort(water.z))
        }
        nz <- length(water.z)
        if (missing(water.col))
            water.col <- oce.colors.gebco(nz, "water", "line")
        if (missing(water.lty))
            water.lty <- rep(par("lty"), nz)
        else if (length(water.lty) == 1)
            water.lty <- rep(water.lty, nz)
        if (missing(water.lwd))
            water.lwd <- rep(par("lwd"), nz)
        else if (length(water.lwd) == 1)
            water.lwd <- rep(water.lwd, nz)
        legend <- c(legend, water.z)
        lwd    <- c(lwd,    water.lwd)
        lty    <- c(lty,    water.lty)
        col    <- c(col,    water.col)
        contour(x$data$longitude, x$data$latitude, x$data$z,
                levels=water.z, lwd=water.lwd, lty=water.lty, col=water.col,
                drawlabels=FALSE, add=TRUE, ...)
    }
    if (zr[2] > 0) {
        if (missing(land.z)) {
            if (zr[1] < 0) {
                land.z <- pretty(c(0, zr[2]))
                land.z <- land.z[land.z!=0]
            } else {
                land.z <- pretty(zr)
            }
        }
        nz <- length(land.z)
        if (nz > 0) {
            if (missing(land.col))
                land.col <- oce.colors.gebco(nz, "land", "line")
            if (missing(land.lty))
                land.lty <- rep(par("lty"), nz)
            else if (length(land.lty) == 1)
                land.lty <- rep(land.lty, nz)
            if (missing(land.lwd))
                land.lwd <- rep(par("lwd"), nz)
            else if (length(land.lwd) == 1)
                land.lwd <- rep(land.lwd, nz)
            legend <- c(legend, land.z)
            lwd    <- c(lwd,    land.lwd)
            lty    <- c(lty,    land.lty)
            col    <- c(col,    land.col)
            contour(x$data$longitude, x$data$latitude, x$data$z,
                    levels=land.z, lwd=land.lwd, lty=land.lty, col=land.col,
                    drawlabels=FALSE, add=TRUE, ...)
        }
    }
    if (!is.null(legend.loc)) {
        o <- rev(order(legend))
        legend(legend.loc, lwd=lwd[o], lty=lty[o],
               bg="white", legend=legend[o], col=col[o])
    }
}

read.topo <- function(file, log.action)
{
    nh <- 6
    header <- readLines(file, n=nh)
    ncols <- as.numeric(strsplit(header[1],"[ ]+",perl=TRUE)[[1]][2])
    nrows <- as.numeric(strsplit(header[2],"[ ]+",perl=TRUE)[[1]][2])
    lon.ll <- as.numeric(strsplit(header[3],"[ ]+",perl=TRUE)[[1]][2])
    lat.ll <- as.numeric(strsplit(header[4],"[ ]+",perl=TRUE)[[1]][2])
    cellsize <- as.numeric(strsplit(header[5],"[ ]+",perl=TRUE)[[1]][2])
    zz <- as.matrix(read.table(file, header=FALSE, skip=nh),byrow=TRUE)
    longitude <- lon.ll + cellsize * seq(0, ncols-1)
    latitude <- lat.ll + cellsize * seq(0, nrows-1)
    z <- t(zz[dim(zz)[1]:1,])
    if (missing(log.action)) log.action <- paste(deparse(match.call()), sep="", collapse="")
    log.item <- processing.log.item(log.action)
    as.topo(longitude, latitude, z, filename=file, log.action=log.item)
}

as.topo <- function(longitude, latitude, z, filename="", log.action)
{
    ncols <- length(longitude)
    nrows <- length(latitude)
    lon.ll <- min(longitude, na.rm=TRUE)
    lat.ll <- min(latitude, na.rm=TRUE)
    dim <- dim(z)
    if (dim[1] != ncols) stop("longitude vector has length ", ncols, ", which does not match matrix width ", dim[1])
    if (dim[2] != nrows) stop("latitude vector has length ", ncols, ", which does not match matrix height ", dim[2])
    data <- list(longitude=longitude, latitude=latitude, z=z)
    metadata <- list(filename=file, ncols=ncols, nrows=nrows, lon.ll=lon.ll, lat.ll=lat.ll)
    if (missing(log.action)) log.action <- paste(deparse(match.call()), sep="", collapse="")
    log.item <- processing.log.item(log.action)
    rval <- list(data=data, metadata=metadata, processing.log=log.item)
    class(rval) <- c("topo", "oce")
    rval
}

summary.topo <- function(object, ...)
{
    if (!inherits(object, "topo")) stop("method is only for topo objects")
    res <- list(lat.range=range(object$data$lat, na.rm=TRUE),
                lon.range=range(object$data$lon, na.rm=TRUE),
                z.range=range(object$data$z, na.rm=TRUE),
                processing.log=processing.log.summary(object))
    class(res) <- "summary.topo"
    res
}

print.summary.topo <- function(x, digits=max(6, getOption("digits") - 1), ...)
{
    cat("\nETOPO dataset\n")
    cat("latitude range:", format(x$lat.range[1], digits),
        " to ", format(x$lat.range[2], digits), "\n")
    cat("longitude range:", format(x$lon.range[1], digits),
        " to ", format(x$lon.range[2], digits), "\n")
    cat("elevation range:", format(x$z.range[1], digits=digits),
        " to ", format(x$z.range[2], digits), "\n")
    print(x$processing.log)
    invisible(x)
}
